####TABLE SEVEN
pkgs <- c("mediation", "dplyr", 
          "haven", "xtable")
sapply(pkgs, require, character.only = T); rm(pkgs)
data <- read_dta("final_dataset.dta")


data$occupation <- 0
data$occupation[data$P2_9 == 99|data$P2_9 == 98] <- NA

data$occupation[data$P2_9 == 4|data$P2_9 == 5] <- 1
data$occupation[data$P2_9 == 6|data$P2_9 == 7|data$P2_9 == 8|data$P2_9 == 9] <- 2
data$occupation[data$P2_9 == 11] <- 3
data$occupation[data$P2_9 == 12] <- 4
data$occupation[data$P2_9 == 13] <- 5
data$occupation[data$P2_9 == 14] <- 6

data$cat_age <- NA
data$cat_age[data$age_at_arrest <= 25] <- 1
data$cat_age[data$age_at_arrest <= 35&data$age_at_arrest > 25] <- 2
data$cat_age[data$age_at_arrest <= 45&data$age_at_arrest > 35] <- 3
data$cat_age[data$age_at_arrest <= 55&data$age_at_arrest > 45] <- 4
data$cat_age[data$age_at_arrest <= 65&data$age_at_arrest > 55] <- 5
data$cat_age[data$age_at_arrest > 65] <- 6

data$education <- NA
data$education[data$P1_24_N == 0] <- 0
data$education[data$P1_24_N == 1|data$P1_24_N == 2] <- 1
data$education[data$P1_24_N == 3|data$P1_24_N == 4|data$P1_24_N == 5] <- 2
data$education[data$P1_24_N == 6|data$P1_24_N == 7] <- 3
data$education[data$P1_24_N == 8|data$P1_24_N == 9] <- 4

data$illiterate <- ifelse(data$P1_22 ==2 | data$P1_23==2, 1, 0)
data$indigenous <- ifelse(data$P1_20 == 1, 1, 0)

data$wealth_q <- NA
data$wealth_q[data$index_wealth < quantile(data$index_wealth, prob = .25, na.rm = T)] <- "0% - 25%"
data$wealth_q[data$index_wealth >= quantile(data$index_wealth, prob = .25, na.rm = T) & data$index_wealth < quantile(data$index_wealth, prob = .5, na.rm = T)] <- "25% - 50%"
data$wealth_q[data$index_wealth >= quantile(data$index_wealth, prob = .5, na.rm = T) & data$index_wealth < quantile(data$index_wealth, prob = .75, na.rm = T)] <- "50% - 75%"
data$wealth_q[data$index_wealth >= quantile(data$index_wealth, prob = .75, na.rm = T)] <- "75%+"

data$muni_pol <- ifelse(
  data$authority2 == "Municipal Police", 1, 0
)
data$min_pol <- ifelse(
  data$authority2 == "Ministerial Police", 1, 0
)

# INST --------------------------------------------------------------------

data_model <- 
  data %>% 
  filter(
    !is.na(cluster)&
      !is.na(year_arrest)&
      !is.na(theft)&
      !is.na(institutional_torture_dum)&
      !is.na(allreform_with_fed)
  ) %>% 
  mutate(
    year_arrest = factor(year_arrest),
    cluster = factor(cluster)
  )

model.mediator <- 
  lm(
    theft ~ allreform_with_fed + year_arrest + cluster,
    data = data_model
  )

model.treatment_1 <- 
  lm(
    institutional_torture_dum ~ theft + allreform_with_fed + year_arrest + cluster,
    data = data_model
  )
cluster_vec <- data_model$cluster
set.seed(2282020)
m_analysis_inst <- 
  mediate(
    model.mediator,
    model.treatment_1,
    treat = "allreform_with_fed",
    mediator = "theft",
    cluster = cluster_vec
  )

# BRUTE FORCE -------------------------------------------------------------

data_model <- 
  data %>% 
  filter(
    !is.na(cluster)&
      !is.na(year_arrest)&
      !is.na(theft)&
      !is.na(brute_force_dum)&
      !is.na(allreform_with_fed)
  ) %>% 
  mutate(
    year_arrest = factor(year_arrest),
    cluster = factor(cluster)
  ) %>% data.frame()
cluster_vec <- data_model$cluster
covar.data <- data_model[,c("year_arrest", "cluster")] %>% data.frame()

model.mediator <- 
  lm(
    theft ~ allreform_with_fed + year_arrest + cluster,
    data = data_model
  )
model.treatment_1 <- 
  lm(
    brute_force_dum ~ 
      theft + allreform_with_fed + year_arrest + cluster,
    data = data_model
  )

set.seed(2292020)
m_analysis_brute <- 
  mediate(
    model.mediator,
    model.treatment_1,
    treat = "allreform_with_fed",
    mediator = "theft",
    cluster = cluster_vec
  )

# THREATS -----------------------------------------------------------------

data_model <- 
  data %>% 
  filter(
    !is.na(data$cluster)&!is.na(data$year_arrest)&
      !is.na(data$theft)&!is.na(data$threat_dum)&!is.na(data$allreform_with_fed)
  ) %>% 
  mutate(
    year_arrest = factor(year_arrest),
    cluster = factor(cluster)
  )

model.mediator <- 
  lm(
    theft ~ allreform_with_fed + year_arrest + cluster,
    data = data_model
  )
model.treatment_1 <- 
  lm(
    threat_dum ~ 
      theft + allreform_with_fed + year_arrest + cluster,
    data = data_model
  )

cluster_vec <- data_model$cluster

covar.data <- data_model[,c("year_arrest", "cluster")] %>% data.frame()

set.seed(3012020)
m_analysis_threat <- 
  mediate(
    model.mediator,
    model.treatment_1,
    treat = "allreform_with_fed",
    mediator = "theft",
    cluster = cluster_vec
  )

theft_models <- 
  list(
    m_analysis_inst
    , m_analysis_brute
    , m_analysis_threat
  )



# MUNPOL ------------------------------------------------------------------

# INST --------------------------------------------------------------------


data_model <- 
  data %>% 
  filter(
    !is.na(data$cluster)&!is.na(data$year_arrest)&!is.na(data$muni_pol)&!is.na(data$institutional_torture_dum)&!is.na(data$allreform_with_fed)
  ) %>% 
  mutate(
    year_arrest = factor(year_arrest),
    cluster = factor(cluster)
  )

model.mediator <- 
  lm(
    muni_pol ~ allreform_with_fed + year_arrest + cluster,
    data = data_model
  )

summary(model.mediator)

model.treatment_1 <- 
  lm(
    institutional_torture_dum ~ muni_pol + allreform_with_fed + year_arrest + cluster,
    data = data_model
  )


cluster_vec <- data_model$cluster

covar.data <- data_model[,c("year_arrest", "cluster")] %>% data.frame()

set.seed(3022020)
m_analysis_inst <- 
  mediate(
    model.mediator,
    model.treatment_1,
    treat = "allreform_with_fed",
    mediator = "muni_pol",
    cluster = cluster_vec
  )
summary(m_analysis_inst)




# BRUTE FORCE -------------------------------------------------------------






data_model <- 
  data %>% 
  filter(
    !is.na(cluster)&
      !is.na(year_arrest)&
      !is.na(muni_pol)&
      !is.na(brute_force_dum)&
      !is.na(allreform_with_fed)
  ) %>% 
  mutate(
    year_arrest = factor(year_arrest),
    cluster = factor(cluster)
  ) %>% data.frame()
cluster_vec <- data_model$cluster
covar.data <- data_model[,c("year_arrest", "cluster")] %>% data.frame()

model.mediator <- 
  lm(
    muni_pol ~ allreform_with_fed + year_arrest + cluster,
    data = data_model
  )
model.treatment_1 <- 
  lm(
    brute_force_dum ~ 
      muni_pol + allreform_with_fed + year_arrest + cluster,
    data = data_model
  )
set.seed(3032020)
m_analysis_brute <- 
  mediate(
    model.mediator,
    model.treatment_1,
    treat = "allreform_with_fed",
    mediator = "muni_pol",
    cluster = cluster_vec
  )

summary(m_analysis_brute)

# THREATS -----------------------------------------------------------------



data_model <- 
  data %>% 
  filter(
    !is.na(cluster)&
      !is.na(year_arrest)&
      !is.na(muni_pol)&
      !is.na(threat_dum)&
      !is.na(allreform_with_fed)
  ) %>% 
  mutate(
    year_arrest = factor(year_arrest),
    cluster = factor(cluster)
  )

model.mediator <- 
  lm(
    muni_pol ~ allreform_with_fed + year_arrest + cluster,
    data = data_model
  )
model.treatment_1 <- 
  lm(
    threat_dum ~ 
      muni_pol + allreform_with_fed + year_arrest + cluster,
    data = data_model
  )

cluster_vec <- data_model$cluster
covar.data <- data_model[,c("year_arrest", "cluster")] %>% data.frame()
set.seed(3042020)

m_analysis_threat <- 
  mediate(
    model.mediator,
    model.treatment_1,
    treat = "allreform_with_fed",
    mediator = "muni_pol",
    cluster = cluster_vec
  )
summary(m_analysis_threat)
police_model <- 
  list(
    m_analysis_inst
    , m_analysis_brute
    , m_analysis_threat
  )

##Torture models
##Police
r_1 <- round(
  c(summary(police_model[[1]])[['d.avg']],
    summary(police_model[[1]])[['d.avg.ci']][1],
    summary(police_model[[1]])[['d.avg.ci']][2]),
  4)
r_2 <- round(
  c(summary(police_model[[1]])[['z.avg']], 
    summary(police_model[[1]])[['z.avg.ci']][1],
    summary(police_model[[1]])[['z.avg.ci']][2]),
  4)
##Theft
r_3 <- round(
  c(summary(theft_models[[1]])[['d.avg']], 
    summary(theft_models[[1]])[['d.avg.ci']][1],
    summary(theft_models[[1]])[['d.avg.ci']][2]),
  4)
r_4 <- round(
  c(summary(theft_models[[1]])[['z.avg']], 
    summary(theft_models[[1]])[['z.avg.ci']][1],
    summary(theft_models[[1]])[['z.avg.ci']][2]),
  4)

###Brute force models
##Police
r_5 <- round(
  c(summary(police_model[[2]])[['d.avg']],
    summary(police_model[[2]])[['d.avg.ci']][1],
    summary(police_model[[2]])[['d.avg.ci']][2]),
  4)
r_6 <- round(
  c(summary(police_model[[2]])[['z.avg']], 
    summary(police_model[[2]])[['z.avg.ci']][1],
    summary(police_model[[2]])[['z.avg.ci']][2]),
  4)
##Theft
r_7 <- round(
  c(summary(theft_models[[2]])[['d.avg']], 
    summary(theft_models[[2]])[['d.avg.ci']][1],
    summary(theft_models[[2]])[['d.avg.ci']][2]),
  4)
r_8 <- round(
  c(summary(theft_models[[2]])[['z.avg']], 
    summary(theft_models[[2]])[['z.avg.ci']][1],
    summary(theft_models[[2]])[['z.avg.ci']][2]),
  4)

###Threat models
##Police
r_9 <- round(
  c(summary(police_model[[3]])[['d.avg']],
    summary(police_model[[3]])[['d.avg.ci']][1],
    summary(police_model[[3]])[['d.avg.ci']][2]),
  4)
r_10 <- round(
  c(summary(police_model[[3]])[['z.avg']], 
    summary(police_model[[3]])[['z.avg.ci']][1],
    summary(police_model[[3]])[['z.avg.ci']][2]),
  4)
##Theft
r_11 <- round(
  c(summary(theft_models[[3]])[['d.avg']], 
    summary(theft_models[[3]])[['d.avg.ci']][1],
    summary(theft_models[[3]])[['d.avg.ci']][2]),
  4)
r_12 <- round(
  c(summary(theft_models[[3]])[['z.avg']], 
    summary(theft_models[[3]])[['z.avg.ci']][1],
    summary(theft_models[[3]])[['z.avg.ci']][2]),
  4)

output <- 
  rbind(
    r_1, r_2,
    r_3, r_4,
    r_5, r_6,
    r_7, r_8,
    r_9, r_10,
    r_11, r_12
  )
colnames(output) <- c(
  "Coefficient",
  "Lower 95 CI",
  "Upper 95 CI"
)

rownames(output) <- 
  c(
    paste(rep(c("Mediator", "Direct"), 2), "torture", sep = ", "),
    paste(rep(c("Mediator", "Direct"), 2), "brute", sep = ", "),
    paste(rep(c("Mediator", "Direct"), 2), "threat", sep = ", ")
  )

xtable(output, digits = 4)
###Significance codes extracted manually using summary() command






